{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 诸闭壳层量子化学方法的密度矩阵"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> 创建时间：2021-01-04；最后修改：2021-06-10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在这份简短笔记中，我们会回顾一些量子化学方法的密度矩阵，及其性质。大体的结论在下述表格中。\n",
    "\n",
    "我们在这里只讨论闭壳层与实函数的情况。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "| 方法 | RHF 轨道基函数 | 能量关系 | $P_p^q$ 对称性 | $\\Gamma_{pr}^{qs}$ 对称性 | 1-RDM 与电子数 | $\\Gamma_{pr}^{qs}$ 与 $P_p^q$ 的关系 | $\\mathbf{P}$ 幂等性 | $P_i^a$ 为零 | $\\Gamma_{ij}^{ab}$ 为零 | $F_p^q$ 对称性 | 1-RDM 偶极矩 |\n",
    "| ------- |:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|\n",
    "| RHF     | √ | √ | √ | √ | √ | √ | √ | √ | √ | √ | √ |\n",
    "| Full-CI | √ | √ | √ | √ | √ | √ | × | × | × | √ | √ |\n",
    "| MP2     | √ | √ | √ | √ | √ | × | × | √ | × | × | × |\n",
    "| CCSD    | √ | √ | √ | √ | √ | √ | × | × | × | × | × |\n",
    "| CISD    | √ | √ | √ | √ | √ | √ | × | × | × | × | × |\n",
    "| CASCI   | × | √ | √ | √ | √ | √ | × | N/A | N/A | × | × |\n",
    "| CASSCF  | × | √ | √ | √ | √ | √ | × | N/A | N/A | √ | √ |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "之所以上面表格中 CASCI、CASSCF 方法不能说 $P_i^a$ (密度矩阵的占据-非占) 与 $\\Gamma_{ij}^{ab}$ (2-RDM 的占据-非占)，是因为它们并是非基于 RHF 参考态的方法，不存在确切的占据与未占轨道。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 预定义"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyscf import gto, scf, mp, cc, ci, mcscf, fci\n",
    "import numpy as np\n",
    "from functools import partial\n",
    "\n",
    "np.einsum = partial(np.einsum, optimize=True)\n",
    "np.set_printoptions(precision=5, linewidth=150, suppress=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "这里讨论的密度矩阵并非是 Full-CI 的情形矩阵，而是会因方法各异而不同的。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "这里采用相对比较严格的 Einstein Summation Convention，即被求和角标必须是一个在上，一个在下。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "这份文档中使用下述上下标：\n",
    "\n",
    "- $p, q, r, s, m, n$ 分子轨道\n",
    "\n",
    "- $i, j$ 分子占据轨道，$a, b$ 分子未占轨道\n",
    "\n",
    "- $\\mu, \\nu, \\kappa, \\lambda$ 原子轨道"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "分子轨道函数 $\\phi_p (\\boldsymbol{r})$ 与原子轨道函数 $\\phi_\\mu (\\boldsymbol{r})$ 之间满足关系 ($C_p^\\mu$ 称为原子轨道系数)\n",
    "\n",
    "$$\n",
    "\\phi_p (\\boldsymbol{r}) = C_p^\\mu \\phi_\\mu (\\boldsymbol{r})\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们假定了研究体系必然是实函数，但我们暂且定义函数的共轭记号如下：(不使用 Einstein Summation)\n",
    "\n",
    "$$\n",
    "\\phi^p (\\boldsymbol{r}) = \\phi_p^* (\\boldsymbol{r})\n",
    "$$\n",
    "\n",
    "分子轨道之间是正交归一的，但原子轨道需要用重叠积分：\n",
    "\n",
    "$$\n",
    "\\int \\phi_p (\\boldsymbol{r}) \\phi^q (\\boldsymbol{r}) \\, \\mathrm{d} \\boldsymbol{r} = \\delta_p^q, \\; \\int \\phi_\\mu (\\boldsymbol{r}) \\phi^\\nu (\\boldsymbol{r}) \\, \\mathrm{d} \\boldsymbol{r} = S_\\mu^\\nu\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Full-CI 密度矩阵的定义与性质"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "这里只作理论上的讨论。Full-CI 密度矩阵程序上的实现会在后面呈现。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "密度矩阵与约化密度的定义有关。由于我们只讨论闭壳层情形，因此波函数可以安全地写成空间坐标的函数。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "### 1-RDM 与基转换关系"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "我们回顾一阶约化密度 $\\rho(\\boldsymbol{r}; \\boldsymbol{r}')$：\n",
    "\n",
    "$$\n",
    "\\rho(\\boldsymbol{r}; \\boldsymbol{r}') = \\idotsint \\Psi^* (\\boldsymbol{r}, \\boldsymbol{r}_2, \\boldsymbol{r}_3, \\cdots, \\boldsymbol{r}_{n_\\mathrm{elec}}) \\Psi (\\boldsymbol{r}', \\boldsymbol{r}_2, \\boldsymbol{r}_3, \\cdots, \\boldsymbol{r}_{n_\\mathrm{elec}}) \\, \\mathrm{d} \\boldsymbol{r}_2 \\, \\mathrm{d} \\boldsymbol{r}_3 \\cdots \\, \\mathrm{d} \\boldsymbol{r}_{n_\\mathrm{elec}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "但是现在只有有限的基函数展开一阶约化密度；如果这组基函数是 RHF 分子轨道 $\\{ \\phi_{p} (\\boldsymbol{r}) \\}$，那么定义下述分子轨道基一阶约化密度矩阵 $P_p^q$ (One-Order Reduced Density Matrix, 1-RDM)\n",
    "\n",
    "$$\n",
    "\\rho(\\boldsymbol{r}; \\boldsymbol{r}') = P_p^q \\phi^p (\\boldsymbol{r}) \\phi_q (\\boldsymbol{r}')\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "如果是原子轨道 $\\{ \\phi_\\mu (\\boldsymbol{r}) \\}$，那么它称为原子轨道基 1-RDM $P_\\mu^\\nu$\n",
    "\n",
    "$$\n",
    "\\rho(\\boldsymbol{r}; \\boldsymbol{r}') = P_\\mu^\\nu \\phi^\\mu (\\boldsymbol{r}) \\phi_\\nu (\\boldsymbol{r}')\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "依据分子轨道与原子轨道间的关系，有\n",
    "\n",
    "$$\n",
    "\\rho(\\boldsymbol{r}; \\boldsymbol{r}') = C_\\mu^p P_p^q C_q^\\nu \\phi^\\mu (\\boldsymbol{r}) \\phi_\\nu (\\boldsymbol{r}')\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "因此，原子轨道基与分子轨道基的 1-RDM 间存在关系\n",
    "\n",
    "$$\n",
    "P_\\mu^\\nu = C_\\mu^p P_p^q C_q^\\nu \n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "### 1-RDM 迹"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "当 $\\boldsymbol{r}, \\boldsymbol{r}'$ 相同时，我们会将一阶约化密度简记为电子密度 $\\rho(\\boldsymbol{r}) = \\rho(\\boldsymbol{r}; \\boldsymbol{r})$。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "一阶约化密度具有积分为电子数的性质：\n",
    "\n",
    "$$\n",
    "\\int \\rho(\\boldsymbol{r}) \\, \\mathrm{d} \\boldsymbol{r} = n_\\mathrm{elec}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "在分子轨道基的表示下，上式可以写为\n",
    "\n",
    "$$\n",
    "P_p^q \\int \\phi^p (\\boldsymbol{r}) \\phi_q (\\boldsymbol{r}) \\, \\mathrm{d} \\boldsymbol{r} = P_p^q \\delta^p_q = \\mathrm{tr} (\\mathbf{P}) = n_\\mathrm{nelec}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "### 1-RDM 对称性"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "首先我们可以证明下述等式：\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "&\\quad\\ \\iint \\phi_p (\\boldsymbol{r}) \\rho(\\boldsymbol{r}; \\boldsymbol{r}') \\phi^q (\\boldsymbol{r}') \\, \\mathrm{d} \\boldsymbol{r} \\, \\mathrm{d} \\boldsymbol{r}' \\\\\n",
    "&= \\iint \\phi_p (\\boldsymbol{r}) P_r^s \\phi^r (\\boldsymbol{r}) \\phi_s (\\boldsymbol{r}') \\phi^q (\\boldsymbol{r}') \\, \\mathrm{d} \\boldsymbol{r} \\, \\mathrm{d} \\boldsymbol{r}' \\\\\n",
    "&= P_r^s \\delta_p^r \\delta_s^q = P_p^q\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "对于上式，如果我们交换被积元变量 $\\boldsymbol{r}, \\boldsymbol{r'}$、并对表达式取共轭，得到 (不使用 Einstein Summation)\n",
    "\n",
    "$$\n",
    "\\iint \\phi_p (\\boldsymbol{r}) \\rho(\\boldsymbol{r}; \\boldsymbol{r}') \\phi^q (\\boldsymbol{r}') \\, \\mathrm{d} \\boldsymbol{r} \\, \\mathrm{d} \\boldsymbol{r}' = \\iint \\phi_q (\\boldsymbol{r}) \\rho^*(\\boldsymbol{r}'; \\boldsymbol{r}) \\phi^p (\\boldsymbol{r}') \\, \\mathrm{d} \\boldsymbol{r} \\, \\mathrm{d} \\boldsymbol{r}'\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "如果我们再利用实数情形下，根据一阶约化密度的定义，有 $\\rho(\\boldsymbol{r}; \\boldsymbol{r}') = \\rho(\\boldsymbol{r}'; \\boldsymbol{r}) = \\rho^*(\\boldsymbol{r}'; \\boldsymbol{r})$，那么可以立即得到 (不使用 Einstein Summation)\n",
    "\n",
    "$$\n",
    "P_p^q = \\iint \\phi_p (\\boldsymbol{r}) \\rho(\\boldsymbol{r}; \\boldsymbol{r}') \\phi^q (\\boldsymbol{r}') \\, \\mathrm{d} \\boldsymbol{r} \\, \\mathrm{d} \\boldsymbol{r}' = \\iint \\phi_q (\\boldsymbol{r}) \\rho(\\boldsymbol{r}; \\boldsymbol{r}') \\phi^p (\\boldsymbol{r}') \\, \\mathrm{d} \\boldsymbol{r} \\, \\mathrm{d} \\boldsymbol{r}' = P_q^p\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "即 1-RDM 矩阵 $P_p^q$ 是对称矩阵。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "### 2-RDM "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "与 1-RDM 相同地，依据二阶约化密度 $\\gamma (\\boldsymbol{r}_1, \\boldsymbol{r}_2; \\boldsymbol{r}'_1, \\boldsymbol{r}'_2)$ 的定义：\n",
    "\n",
    "$$\n",
    "\\gamma (\\boldsymbol{r}_1, \\boldsymbol{r}_2; \\boldsymbol{r}'_1, \\boldsymbol{r}'_2) = \\idotsint \\Psi^* (\\boldsymbol{r}_1, \\boldsymbol{r}_2, \\boldsymbol{r}_3, \\cdots, \\boldsymbol{r}_{n_\\mathrm{elec}}) \\Psi (\\boldsymbol{r}'_1, \\boldsymbol{r}'_2, \\boldsymbol{r}_3, \\cdots, \\boldsymbol{r}_{n_\\mathrm{elec}}) \\, \\mathrm{d} \\boldsymbol{r}_3 \\cdots \\, \\mathrm{d} \\boldsymbol{r}_{n_\\mathrm{elec}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "用分子轨道基作展开，可以定义分子轨道基的二阶约化密度矩阵 $\\Gamma_{pr}^{qs}$ (Two-Order Reduced Density Matrix, 2-RDM)\n",
    "\n",
    "$$\n",
    "\\gamma (\\boldsymbol{r}_1, \\boldsymbol{r}_2; \\boldsymbol{r}'_1, \\boldsymbol{r}'_2) = \\Gamma_{pr}^{qs} \\phi^p (\\boldsymbol{r}_1) \\phi_q (\\boldsymbol{r}'_1) \\phi^r (\\boldsymbol{r}_2) \\phi_s (\\boldsymbol{r}'_2)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "原子与分子轨道基转换也与 1-RDM 类似：\n",
    "\n",
    "$$\n",
    "\\Gamma_{\\mu \\kappa}^{\\nu \\lambda} = C_\\mu^p C^\\nu_q \\Gamma_{pr}^{qs} C_\\kappa^r C^\\lambda_s\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "### 2-RDM 与 1-RDM 的关系"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "出于全同粒子的性质，2-RDM 与 1-RDM 之间存在关系：\n",
    "\n",
    "$$\n",
    "\\rho(\\boldsymbol{r}_1; \\boldsymbol{r}'_1) = \\frac{1}{n_\\mathrm{elec} - 1} \\iint \\gamma (\\boldsymbol{r}_1, \\boldsymbol{r}_2; \\boldsymbol{r}'_1, \\boldsymbol{r}_2) \\, \\mathrm{d} \\boldsymbol{r}_2\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "对上式展开并作一部分积分后，可以得到\n",
    "\n",
    "$$\n",
    "P_p^q \\phi^p (\\boldsymbol{r}_1) \\phi_q (\\boldsymbol{r}'_1) = \\frac{1}{n_\\mathrm{elec} - 1} \\Gamma_{pr}^{qm} \\phi^p (\\boldsymbol{r}_1) \\phi_q (\\boldsymbol{r}'_1) \\delta_m^r\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "由于上式要在任意的 $\\boldsymbol{r}_1, \\boldsymbol{r}'_1$ 的取值下成立，因此可以认为\n",
    "\n",
    "$$\n",
    "P_p^q = \\frac{1}{n_\\mathrm{elec} - 1} \\Gamma_{pr}^{qm} \\delta_m^r\n",
    "$$\n",
    "\n",
    "注意上式要对等式右边作关于 $r, m$ 角标的求和。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "### 2-RDM 对称性"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "分析 2-RDM 对称性相对比较麻烦。这里就略过讨论了。我们仅指出，实数闭壳层下的 2-RDM 应当具有二重对称性，而不具有更高的对称性：\n",
    "\n",
    "$$\n",
    "\\Gamma_{pr}^{qs} = \\Gamma_{rp}^{sq}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "### 密度矩阵与能量"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "这是最关键的一个性质。密度矩阵可以用来表示电子态的能量。\n",
    "\n",
    "现在记原子轨道基组下的单电子算符积分为 $h_\\mu^\\nu$、双电子算符积分为 $g_{\\mu \\kappa}^{\\nu \\lambda}$，其中单电子算符包含动能、原子核-电子库伦势能、电场势能等贡献，双电子算符包含电子-电子库伦势能贡献。那么，体系单点能为\n",
    "\n",
    "$$\n",
    "E_\\mathrm{tot} = E_\\mathrm{elec} + E_\\mathrm{nuc} = h_\\mu^\\nu P_\\nu^\\mu + \\frac{1}{2} g_{\\mu \\kappa}^{\\nu \\lambda} \\Gamma_{\\nu \\lambda}^{\\mu \\kappa} + E_\\mathrm{nuc}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "### 1-RDM 与偶极矩"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Full-CI 的 1-RDM 可以直接用以计算偶极矩；以 $z$ 轴方向施加的电场为例：\n",
    "\n",
    "$$\n",
    "d_z = - z_\\mu^\\nu P_\\nu^\\mu\n",
    "$$\n",
    "\n",
    "其中 $z_\\mu^\\nu = \\langle \\mu | z | \\nu \\rangle$。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 广义 Fock 矩阵"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "广义 Fock 矩阵定义为\n",
    "\n",
    "$$\n",
    "F_p^q = h_p^r P_r^q + g_{pr}^{ms} \\Gamma_{ms}^{qr}\n",
    "$$\n",
    "\n",
    "特别地，在 RHF 下，它是对角矩阵。而在 Full-CI 下，它会是对称矩阵。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## RHF 密度矩阵特有性质"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 幂等性"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在 Conanical-HF 下，幂等性是几乎显然的：$P_p^q$ 一定是对角矩阵；如果 $p$ 所代表的轨道是占据轨道，那么一定填了因为填了两个电子而值为 2，否则为零。因此，Conanical-HF 下一定满足\n",
    "\n",
    "$$\n",
    "P_p^m P_m^q = 2 P_p^q\n",
    "$$\n",
    "\n",
    "一般程序都只会给出 Conanical-HF 的结果。但若讨论 Nonconanical-HF 时 $P_p^q$ 未必是对角矩阵，但上述结论应仍然成立。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 非占-占据部分为零"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hartree-Fock 方法严格地将轨道分为占据与非占据。因此，Canonical 或 Nonconanical HF 方法都会保证 1-RDM 是块状对角化的；即在占据-非占 $P_i^a$，非占-占据 $P_a^i$，非占-非占 $P_a^b$ 均严格为零。对于 2-RDM 是类似的。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "由于 Hartree-Fock 方法没有考虑非占轨道的贡献，因此任何 Post-HF 方法均一定程度上有激发态的贡献。一般来说，非占-非占的 $P_a^b$ 贡献总是存在的；但占据-非占或非占-占据的 $P_i^a$ 与 $P_a^i$ 则未必存在。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 通用计算函数"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "下面的文档仅仅是验证开头的表格的代码。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 分子定义：水分子"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- `mol` 水分子实例；\n",
    "\n",
    "- `nelec` $n_\\mathrm{elec}$ 电子数；\n",
    "\n",
    "- `nocc` $n_\\mathrm{occ}$ 占据轨道数；\n",
    "\n",
    "- `h` 原子轨道基 $h_\\mu^\\nu$，维度 $(\\mu, \\nu)$；\n",
    "\n",
    "- `g` 原子轨道基 $g_{\\mu \\kappa}^{\\nu \\lambda}$，维度 $(\\mu, \\nu, \\kappa, \\lambda)$；\n",
    "\n",
    "- `S` 原子轨道基 $S_\\mu^\\nu$，维度 $(\\mu, \\nu)$；\n",
    "\n",
    "- `mf_rhf` RHF 实例。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<pyscf.gto.mole.Mole at 0x7feae8ff9580>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mol = gto.Mole()\n",
    "mol.atom = \"\"\"\n",
    "O  0. 0. 0.\n",
    "H  0. 0. 1.\n",
    "H  0. 1. 0.\n",
    "\"\"\"\n",
    "mol.basis = \"6-31G\"\n",
    "mol.verbose = 0\n",
    "mol.build()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(10, 5)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nelec = mol.nelectron\n",
    "nocc = mol.nelec[0]\n",
    "nelec, nocc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "h = mol.intor(\"int1e_kin\") + mol.intor(\"int1e_nuc\")\n",
    "g = mol.intor(\"int2e\")\n",
    "S = mol.intor(\"int1e_ovlp\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "mf_rhf = scf.RHF(mol).run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 验证能量表达式"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "验证\n",
    "\n",
    "$$\n",
    "E_\\mathrm{tot} = h_\\mu^\\nu P_\\nu^\\mu + \\frac{1}{2} g_{\\mu \\kappa}^{\\nu \\lambda} \\Gamma_{\\nu \\lambda}^{\\mu \\kappa} + E_\\mathrm{nuc} = h_p^q P_q^p + \\frac{1}{2} g_{pr}^{qs} \\Gamma_{qs}^{pr} + E_\\mathrm{nuc}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "eng_nuc = mol.energy_nuc()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def verify_energy_relation(eng, eng_nuc, rdm1, rdm2, h_mo, g_mo):\n",
    "    return np.allclose(np.einsum(\"pq, qp ->\", h_mo, rdm1) + 0.5 * np.einsum(\"pqrs, qpsr ->\", g_mo, rdm2) + eng_nuc, eng)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 验证 1-RDM 对称性"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "验证 $P_p^q = P_q^p$。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def verify_rdm1_symm(rdm1):\n",
    "    # Output: 1-RDM symmetric property\n",
    "    return np.allclose(rdm1, rdm1.T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 验证 2-RDM 对称性"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "验证 $\\Gamma_{pr}^{qs} = \\Gamma_{rp}^{sq}$。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def verify_rdm2_symm(rdm2):\n",
    "    return np.allclose(rdm2, np.einsum(\"pqrs -> rspq\", rdm2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 验证 1-RDM 的迹"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "验证 $P_p^r \\delta_r^p = n_\\mathrm{elec}$。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def verify_rdm1_tr(rdm1):\n",
    "    return np.allclose(rdm1.trace(), nelec)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "### 验证 1-RDM 与 2-RDM 的关系"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "验证 $P_p^q = (n_\\mathrm{elec} - 1)^{-1} \\Gamma_{pr}^{qm} \\delta_m^r$。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "def verify_rdm12_relation(rdm1, rdm2):\n",
    "    return np.allclose(rdm1, (nelec - 1)**-1 * rdm2.diagonal(axis1=-1, axis2=-2).sum(axis=-1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 验证 1-RDM 幂等性"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "验证 $P_p^m P_m^q = 2 P_p^q$。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def verify_rdm1_idomp(rdm1):\n",
    "    return np.allclose(rdm1 @ rdm1, 2 * rdm1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "### 验证 $P_i^a$ 为零"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "这里实际上同时验证 $P_a^i$ 是否为零。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "def verify_rdm1_ov(rdm1):\n",
    "    mat1 = rdm1[nocc:, :nocc]\n",
    "    mat2 = rdm1[:nocc, nocc:]\n",
    "    return np.allclose(mat1, np.zeros_like(mat1)) and np.allclose(mat2, np.zeros_like(mat2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "### 验证 $\\Gamma_{ij}^{ab}$ 为零"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": [
    "def verify_rdm2_ovov(rdm2):\n",
    "    mat = rdm2[:nocc, nocc:, :nocc, nocc:]\n",
    "    return np.allclose(mat, np.zeros_like(mat))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 验证广义 Fock 矩阵对称性"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "F_p^q = h_p^r P_r^q + g_{pr}^{ms} \\Gamma_{ms}^{qr}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def verify_gF_symm(rdm1, rdm2, h_mo, g_mo):\n",
    "    gF = np.einsum(\"pr, rq -> pq\", h_mo, rdm1) + np.einsum(\"pmrs, mqsr -> pq\", g_mo, rdm2)\n",
    "    return np.allclose(gF, gF.T, atol=1e-4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 偶极矩的验证"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "通过 1-RDM 计算的偶极矩为 (不考虑原子核影响)\n",
    "\n",
    "$$\n",
    "d_z = - z_\\mu^\\nu P_\\nu^\\mu\n",
    "$$\n",
    "\n",
    "但另一种偶极矩的计算方式是对 $h_\\mu^\\nu$ 作更改，求得该情形下的能量作数值差分得到。数值差分的间隙设定为 1e-4 单位电场强度。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "h_field = 1e-4\n",
    "\n",
    "def get_hcore_p(mol_=mol):\n",
    "    return mol.intor(\"int1e_kin\") + mol.intor(\"int1e_nuc\") - h_field * mol.intor(\"int1e_r\")[2]\n",
    "def get_hcore_m(mol_=mol):\n",
    "    return mol.intor(\"int1e_kin\") + mol.intor(\"int1e_nuc\") + h_field * mol.intor(\"int1e_r\")[2]\n",
    "\n",
    "mf_rhf_p, mf_rhf_m = scf.RHF(mol), scf.RHF(mol)\n",
    "mf_rhf_p.get_hcore = get_hcore_p\n",
    "mf_rhf_m.get_hcore = get_hcore_m\n",
    "mf_rhf_p.run(), mf_rhf_m.run()\n",
    "\n",
    "charges = mol.atom_charges()\n",
    "coords  = mol.atom_coords()\n",
    "nucl_dip = np.einsum('i,ix->x', charges, coords)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "def verify_dip(method, rdm1, z_intg):\n",
    "    mf_met_m, _, _, _ = method(mf_rhf_m)\n",
    "    mf_met_p, _, _, _ = method(mf_rhf_p)\n",
    "    dip_num = (mf_met_p.e_tot - mf_met_m.e_tot) / (2 * h_field) + nucl_dip[2]\n",
    "    dip_rdm1 = - (rdm1 * z_intg).sum() + nucl_dip[2]\n",
    "    return np.allclose(dip_num, dip_rdm1, atol=1e-4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 各种方法的验证"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 总验证程序"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "def verify_all(method):\n",
    "    # rdm1, rdm2 here are both in mo_basis\n",
    "    mf_met, C, rdm1, rdm2 = method(mf_rhf)\n",
    "    h_mo = C.T @ h @ C\n",
    "    g_mo = np.einsum(\"up, vq, uvkl, kr, ls -> pqrs\", C, C, g, C, C)\n",
    "    z_intg = C.T @ mol.intor(\"int1e_r\")[2] @ C\n",
    "    print(\"===  Energy Relat  ===  \", verify_energy_relation(mf_met.e_tot, eng_nuc, rdm1, rdm2, h_mo, g_mo))\n",
    "    print(\"===   1-RDM Symm   ===  \", verify_rdm1_symm(rdm1))\n",
    "    print(\"===   2-RDM Symm   ===  \", verify_rdm2_symm(rdm2))\n",
    "    print(\"===   1-RDM Trace  ===  \", verify_rdm1_tr(rdm1))\n",
    "    print(\"===  12-RDM Relat  ===  \", verify_rdm12_relation(rdm1, rdm2))\n",
    "    print(\"===   1-RDM Idomp  ===  \", verify_rdm1_idomp(rdm1))\n",
    "    print(\"===   1-RDM ov     ===  \", verify_rdm1_ov(rdm1))\n",
    "    print(\"===   2-RDM ovov   ===  \", verify_rdm2_ovov(rdm2))\n",
    "    print(\"=== GenFock Symm   ===  \", verify_gF_symm(rdm1, rdm2, h_mo, g_mo))\n",
    "    print(\"===   1-RDM Dipole ===  \", verify_dip(method, rdm1, z_intg))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### RHF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "def method_rhf(mf_rhf):\n",
    "    mf_met = mf_rhf\n",
    "    C = mf_rhf.mo_coeff\n",
    "    Cinv = np.linalg.inv(C)\n",
    "    # In AO basis\n",
    "    rdm1 = mf_rhf.make_rdm1()\n",
    "    rdm2 = np.einsum(\"uv, kl -> uvkl\", rdm1, rdm1) - 0.5 * np.einsum(\"uv, kl -> ukvl\", rdm1, rdm1)\n",
    "    # Transform to MO basis\n",
    "    rdm1 = np.einsum(\"pu, uv, qv -> pq\", Cinv, rdm1, Cinv)\n",
    "    rdm2 = np.einsum(\"pu, qv, uvkl, rk, sl -> pqrs\", Cinv, Cinv, rdm2, Cinv, Cinv)\n",
    "    return mf_met, C, rdm1, rdm2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "===  Energy Relat  ===   True\n",
      "===   1-RDM Symm   ===   True\n",
      "===   2-RDM Symm   ===   True\n",
      "===   1-RDM Trace  ===   True\n",
      "===  12-RDM Relat  ===   True\n",
      "===   1-RDM Idomp  ===   True\n",
      "===   1-RDM ov     ===   True\n",
      "===   2-RDM ovov   ===   True\n",
      "=== GenFock Symm   ===   True\n",
      "===   1-RDM Dipole ===   True\n"
     ]
    }
   ],
   "source": [
    "verify_all(method_rhf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Full-CI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "def method_fci(mf_rhf):\n",
    "    mf_met = fci.FCI(mf_rhf).run()\n",
    "    C = mf_rhf.mo_coeff\n",
    "    # In MO basis\n",
    "    rdm1, rdm2 = mf_met.make_rdm12(mf_met.ci, mol.nao, mol.nelec)\n",
    "    return mf_met, C, rdm1, rdm2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "===  Energy Relat  ===   True\n",
      "===   1-RDM Symm   ===   True\n",
      "===   2-RDM Symm   ===   True\n",
      "===   1-RDM Trace  ===   True\n",
      "===  12-RDM Relat  ===   True\n",
      "===   1-RDM Idomp  ===   False\n",
      "===   1-RDM ov     ===   False\n",
      "===   2-RDM ovov   ===   False\n",
      "=== GenFock Symm   ===   True\n",
      "===   1-RDM Dipole ===   True\n"
     ]
    }
   ],
   "source": [
    "verify_all(method_fci)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### MP2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "def method_mp2(mf_rhf):\n",
    "    mf_met = mp.MP2(mf_rhf).run()\n",
    "    C = mf_rhf.mo_coeff\n",
    "    # In MO basis\n",
    "    rdm1, rdm2 = mf_met.make_rdm1(), mf_met.make_rdm2()\n",
    "    return mf_met, C, rdm1, rdm2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "===  Energy Relat  ===   True\n",
      "===   1-RDM Symm   ===   True\n",
      "===   2-RDM Symm   ===   True\n",
      "===   1-RDM Trace  ===   True\n",
      "===  12-RDM Relat  ===   False\n",
      "===   1-RDM Idomp  ===   False\n",
      "===   1-RDM ov     ===   True\n",
      "===   2-RDM ovov   ===   False\n",
      "=== GenFock Symm   ===   False\n",
      "===   1-RDM Dipole ===   False\n"
     ]
    }
   ],
   "source": [
    "verify_all(method_mp2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### CCSD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "def method_ccsd(mf_rhf):\n",
    "    mf_met = cc.CCSD(mf_rhf).run()\n",
    "    C = mf_rhf.mo_coeff\n",
    "    # In MO basis\n",
    "    rdm1, rdm2 = mf_met.make_rdm1(), mf_met.make_rdm2()\n",
    "    return mf_met, C, rdm1, rdm2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "===  Energy Relat  ===   True\n",
      "===   1-RDM Symm   ===   True\n",
      "===   2-RDM Symm   ===   True\n",
      "===   1-RDM Trace  ===   True\n",
      "===  12-RDM Relat  ===   True\n",
      "===   1-RDM Idomp  ===   False\n",
      "===   1-RDM ov     ===   False\n",
      "===   2-RDM ovov   ===   False\n",
      "=== GenFock Symm   ===   False\n",
      "===   1-RDM Dipole ===   False\n"
     ]
    }
   ],
   "source": [
    "verify_all(method_ccsd)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### CISD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "def method_cisd(mf_rhf):\n",
    "    mf_met = ci.CISD(mf_rhf).run()\n",
    "    C = mf_rhf.mo_coeff\n",
    "    # In MO basis\n",
    "    rdm1, rdm2 = mf_met.make_rdm1(), mf_met.make_rdm2()\n",
    "    return mf_met, C, rdm1, rdm2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "===  Energy Relat  ===   True\n",
      "===   1-RDM Symm   ===   True\n",
      "===   2-RDM Symm   ===   True\n",
      "===   1-RDM Trace  ===   True\n",
      "===  12-RDM Relat  ===   True\n",
      "===   1-RDM Idomp  ===   False\n",
      "===   1-RDM ov     ===   False\n",
      "===   2-RDM ovov   ===   False\n",
      "=== GenFock Symm   ===   False\n",
      "===   1-RDM Dipole ===   False\n"
     ]
    }
   ],
   "source": [
    "verify_all(method_cisd)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### CASCI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "def method_casci(mf_rhf):\n",
    "    mf_met = mcscf.CASCI(mf_rhf, ncas=4, nelecas=4).run()\n",
    "    C = mf_met.mo_coeff\n",
    "    Cinv = np.linalg.inv(C)\n",
    "    # In AO basis\n",
    "    rdm1, rdm2 = mcscf.addons.make_rdm12(mf_met)\n",
    "    # Transform to MO basis\n",
    "    rdm1 = np.einsum(\"pu, uv, qv -> pq\", Cinv, rdm1, Cinv)\n",
    "    rdm2 = np.einsum(\"pu, qv, uvkl, rk, sl -> pqrs\", Cinv, Cinv, rdm2, Cinv, Cinv)\n",
    "    return mf_met, C, rdm1, rdm2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "===  Energy Relat  ===   True\n",
      "===   1-RDM Symm   ===   True\n",
      "===   2-RDM Symm   ===   True\n",
      "===   1-RDM Trace  ===   True\n",
      "===  12-RDM Relat  ===   True\n",
      "===   1-RDM Idomp  ===   False\n",
      "===   1-RDM ov     ===   False\n",
      "===   2-RDM ovov   ===   False\n",
      "=== GenFock Symm   ===   False\n",
      "===   1-RDM Dipole ===   False\n"
     ]
    }
   ],
   "source": [
    "verify_all(method_casci)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### CASSCF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "def method_casscf(mf_rhf):\n",
    "    mf_met = mcscf.CASSCF(mf_rhf, ncas=4, nelecas=4).run()\n",
    "    C = mf_met.mo_coeff\n",
    "    Cinv = np.linalg.inv(C)\n",
    "    # In AO basis\n",
    "    rdm1, rdm2 = mcscf.addons.make_rdm12(mf_met)\n",
    "    # Transform to MO basis\n",
    "    rdm1 = np.einsum(\"pu, uv, qv -> pq\", Cinv, rdm1, Cinv)\n",
    "    rdm2 = np.einsum(\"pu, qv, uvkl, rk, sl -> pqrs\", Cinv, Cinv, rdm2, Cinv, Cinv)\n",
    "    return mf_met, C, rdm1, rdm2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "===  Energy Relat  ===   True\n",
      "===   1-RDM Symm   ===   True\n",
      "===   2-RDM Symm   ===   True\n",
      "===   1-RDM Trace  ===   True\n",
      "===  12-RDM Relat  ===   True\n",
      "===   1-RDM Idomp  ===   False\n",
      "===   1-RDM ov     ===   False\n",
      "===   2-RDM ovov   ===   False\n",
      "=== GenFock Symm   ===   True\n",
      "===   1-RDM Dipole ===   True\n"
     ]
    }
   ],
   "source": [
    "verify_all(method_casscf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 补充"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "感谢 [hebrewsnabla](https://github.com/hebrewsnabla) 对 CASCI、CASSCF 密度矩阵的讨论。"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
